"""
We use the SymPy Python module, available at http://github.com/sympy (offered under BSD license).

We compute the expected total reward on MRMs. For this we use:
    - the Value Iteration Algorithm, [Puterman. Markov Decision Processes: Discrete Stochastic Dynamic Programming, 1994, Chapter 6.3.2] 
        - this does not converge well when we have symbolic transition probabilities.
    - Algorithm 2 from the [Hahn_SPIN2009] paper:
        Ernst Moritz Hahn, Holger Hermanns and Lijun Zhang. "Probabilistic Reachability for Parametric 
            Markov Models", SPIN 2009

!!!!TODO:
   ComputeTotalExpectedRewardWithValueIteration(): we currently use only the state rewards, not transition rewards.
   for method ComputeTotalExpectedRewardWithAlgo2_SPIN2009() - make param sB a list.
**
DONE:
    Make class MDP --> MRM (since we currently work with only 1 action).
    PreXPost
    transition costs besides state costs - rename also such that we have name r (rewards) instead of costs.
    remove a state from the MRM state space
"""


"""
N = 100
#20 iter
#a = 1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N

a = 1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(1.0 + 1.0*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N)*(N - 1)/N

print a
"""

import random
import sys
import traceback

import sympy
#import sympy.core.numbers
#We need to import the symbols like this to check them with isinstance().
from sympy.core.numbers import Zero, Float
from sympy.core.numbers import Mul, Pow, Add
import numpy


MY_DEBUG_STDOUT = True
MY_DEBUG_STDERR = True

#LAMBDA = 0.99 #discount factor
LAMBDA = 1.0 #discount factor - for symbolic computation we can use it
NUM_MAX_ITERATIONS_VALUE_ITERATION = 10 #discount factor

EPS_SMALL = 1.0E-6

class MRM:
    #NUM_STATES_MDP_MAX = 100
    #NUM_ACTIONS_MAX = 10

    numStateVars = None
    cardinality = None
    numReachableStates = None
    #numActions = None

    transitionMatrix = None
    rewardState = None
    rewardTransition = None

    selfPreXPost = None
    pre = None
    post = None

    crtState = None #It is factored representation.

    v = None #value function

    #def PrintTransitionMatrix(self):
    def PrintMatrix(self, matrix):
        #if self.numActions > 1
        #print "PrintMatrix(): printing only for action 0"
        #print "Entered PrintMatrix(matrix = %d):" % (matrix)
        #print "Entered PrintMatrix(matrix = %s):" % (str(hex(id(matrix))))
        print "Entered PrintMatrix(matrix = 0x%x):" % (id(matrix))

        for i in range(self.numReachableStates):
            for j in range(self.numReachableStates):
                if (type(matrix[i][j]) is float):
                    print "%.1f" % matrix[i][j] + "    ",
                elif (isinstance(matrix[i][j], Float)):
                    print "%.1f" % matrix[i][j] + "    ",
                #elif (isinstance(matrix[i][j], sympy.core.numbers.Float)):
                #    print "%.1f" % matrix[i][j] + "    ",
                #elif ((not isinstance(matrix[i][j], sympy.core.mul.Mul)) and 
                #    (not isinstance(matrix[i][j], sympy.core.power.Pow))):
                elif ((not isinstance(matrix[i][j], Mul)) and #AttributeError: 'module' object has no attribute 'mul'
                    (not isinstance(matrix[i][j], Pow)) and
                    (not isinstance(matrix[i][j], Add))
                ):
                    print "%.1f" % matrix[i][j] + "    ",
                #sympy.core.numbers.Zero
                #elif (type(matrix[i][j]) is double):
                #    print "%.1f" % matrix[i][j] + "   !",
                else:
                    print str(matrix[i][j]) + "    ",
                #print type(matrix[i][j])
                #print "PrintTransitionMatrix(): self.transitionMatrix[%d][%d] = %s" % (i, j, str(self.transitionMatrix[i][j]))
            print ""
            sys.stdout.flush()
        print ""


    def __init__(self, numStateVars, stateVarsCardinality, initState):
        self.numStateVars = numStateVars
        self.cardinality = stateVarsCardinality

        #Initial state
        self.crtState = initState

        self.numReachableStates = 1
        for i in range(self.numStateVars):
            self.numReachableStates *= self.cardinality[i]

        #self.numActions = numActions

        #Initializing the value function array
        """
        [Puterman, Chapter 4.1.2]
        Let v_N^{\pi}(s) represent the expected total reward over the decision making horizon if
          policy \pi is used and the system is in state s at the first decision epoch.
        For $\pi in \pi^{HR}$, it is defined by:
        """
        """
        self.v = numpy.array([0]) #, order='C')
        self.v.resize((NUM_MAX_ITERATIONS_VALUE_ITERATION, self.numReachableStates))
        """
        #We can put here also:
        self.v = [None] * NUM_MAX_ITERATIONS_VALUE_ITERATION
        for n in range(NUM_MAX_ITERATIONS_VALUE_ITERATION):
            self.v[n] = [0] * self.numReachableStates
        #print "self.v =", self.v

        #See Example 6.3.1 (page 164) - they also set v[0][*] to 0
        for s in range(self.numReachableStates):
            self.v[0][s] = 0.0

        self.rewardTransition = [None] * self.numReachableStates
        for i in range(self.numReachableStates):
            self.rewardTransition[i] = [0.0] * self.numReachableStates

        #Initializing the transition probability matrix of the MRM (for each action).
        """
        Since we can have symbolic transition probabilities we don't use this:
        self.transitionMatrix = numpy.array([0.0]) #, order='C')
        self.transitionMatrix.resize((self.numActions, self.numReachableStates, self.numReachableStates))
        """

        """
        self.transitionMatrix = [None] * self.numActions
        for i in range(self.numActions):
            self.transitionMatrix[i] = [None] * self.numReachableStates
            for j in range(self.numReachableStates):
                self.transitionMatrix[i][j] = [0.0] * self.numReachableStates
        """
        self.transitionMatrix = [None] * self.numReachableStates
        for i in range(self.numReachableStates):
            self.transitionMatrix[i] = [0.0] * self.numReachableStates


    def ComputeTotalExpectedRewardWithValueIteration(self):
        #print "self.transitionMatrix =", self.transitionMatrix
        self.PrintMatrix(self.transitionMatrix)

        print "self.rewardState =", self.rewardState

        #Value Iteration Algorithm, [Puterman, Chapter 6.3.2]:

        """
        We cannot apply the termination condition as in [Puterman, Chapter 6.3.2], since we deal with symbolic values, as well.
        Therefore, we use as termination condition if the current v[n + 1] and previous one v[n] are equal.
        """
        def Compute_Max_a(s):
            myMax = -1000
            #myMax = 0

            """
            This is basically the implementation of step (6.3.2) from Value Iteration Algorithm, [Puterman, Chapter 6.3.2]:
            """
            sumAll = self.rewardState[s]
            for j in range(self.numReachableStates):
                sumAll += LAMBDA * self.transitionMatrix[s][j] * self.v[n][j]
                print "j = %d, sumAll = %s" % (j, sumAll)

            #myMax = sumAll
            #"""
            if (myMax < sumAll):
                myMax = sumAll
            #"""
            return myMax

        for n in range(NUM_MAX_ITERATIONS_VALUE_ITERATION - 1):
            for s in range(self.numReachableStates):
                res = Compute_Max_a(s)
                res = sympy.simplify(res)

                print "s =", s
                print "n =", n
                print "res =", res
                sys.stdout.flush()
                self.v[n + 1][s] = res
            for s in range(self.numReachableStates):
                if (self.v[n + 1][s] != self.v[n][s]):
                    break
            else:
                #Only executed if break was not used.
                #If we reach here, it means that v[n + 1] == v[n], for all s, and we don't continue anymore.
                lastOne = n
                break
        else:
            lastOne = NUM_MAX_ITERATIONS_VALUE_ITERATION - 1

        print "self.v =", self.v
        print "lastOne =", lastOne
        print "self.v_n =", self.v[lastOne]


    ###############################################################################
    ###############################################################################
    ################ComputeTotalExpectedRewardWithAlgo2_SPIN2009###################
    ###############################################################################
    ###############################################################################

    def ComputePresAndPosts(self):
        """
        #This is not good - this is a serious mistake one can do in Python - when using * the same object reference is used for all instances...
        self.pre = [[]] * self.numReachableStates
        self.post = [[]] * self.numReachableStates
        """
        self.pre = []
        self.post = []
        for s in range(self.numReachableStates):
            self.pre.append([])
            self.post.append([])

        #for s in range(self.numReachableStates):
        #    print "     pre(%d) = %s" % (s, str(self.pre[s]))
        #    print "     post(%d) = %s" % (s, str(self.post[s]))

        for s in range(self.numReachableStates):
            #self.pre[s] = []
            #self.post[s] = []
            for i in range(self.numReachableStates):
                if ( (abs(self.transitionMatrix[s][i] - 0.0) > EPS_SMALL) 
                    and (i != s) ): #For ComputeTotalExpectedRewardWithAlgo2_SPIN2009 to work we don't consider as pre and post the state s itself, if it has a self-loop.
                    #print "PreXPost(): self.transitionMatrix[%d][%d] = %s" % (i, s1, str(self.transitionMatrix[i][s1]))
                    self.post[s].append(i)
                    self.pre[i].append(s)
                    #print "s = %d, i = %d" % (s, i)
                    #print "     pre(%d) = %s" % (s, str(self.pre[s]))
                    #print "     post(%d) = %s" % (i, str(self.post[i]))
            #for i in range(self.numReachableStates):
            #    if (abs(self.transitionMatrix[s][i] - 0.0) > EPS_SMALL):
            #        self.post[i].append(s)

        for s in range(self.numReachableStates):
            print "     pre(%d) = %s" % (s, str(self.pre[s]))
            print "     post(%d) = %s" % (s, str(self.post[s]))


    def PreXPost(self, s1, s2):
        """
        pre_s1 = []
        post_s2 = []

        #self.PrintMatrix(self.transitionMatrix)

        for i in range(self.numReachableStates):
            if (abs(self.transitionMatrix[i][s1] - 0.0) > EPS_SMALL):
                #print "PreXPost(): self.transitionMatrix[%d][%d] = %s" % (i, s1, str(self.transitionMatrix[i][s1]))
                pre_s1.append(i)
            if (abs(self.transitionMatrix[s2][i] - 0.0) > EPS_SMALL):
                post_s2.append(i)
        """

        res = []
        for i in self.pre[s1]:
            for j in self.post[s2]:
                res.append( (i, j) )

        print "PreXPost(%d, %d): pre(%d) = %s" % (s1, s2, s1, str(self.pre[s1]))
        print "                  post(%d) = %s" % (s2, str(self.post[s2]))
        print "                  res = %s" % (str(res))
        return res


    def ComputeSelfPreXPost(self):
        self.selfPreXPost = [None] * self.numReachableStates
        for s in range(self.numReachableStates):
            self.selfPreXPost[s] = self.PreXPost(s, s)


    def ComputeTotalExpectedRewardWithAlgo2_SPIN2009(self, s0, sB):
        #Note that the sB state is absorbing, as the Algo specifies.

        #Applying Algortithm 2 from [SPIN09]
        """
        Example from [SPIN09] paper:
          Initially:
            P(s1, s2) = pd
            P(s1, s) = pa
            P(s, s2) = pb
            P(s, s) = pc
            P(s1, s1) = 0??
            P(s2, s2) = 0??

          We eliminate state s and apply formula from Algo 1:
            P(s1, s2) = P(s1, s2) + P(s1, s) * 1 / (1 - pc) * P(s, s2) = pd + pa * 1/(1-pc) * pb
        """

        #Precompute the PreXPost(s,s), since we will change the transition probability in ComputeTotalExpectedRewardWithAlgo2_SPIN2009().
        self.ComputePresAndPosts()
        self.ComputeSelfPreXPost()

        #self.re = [0] * self.numReachableStates

        #self.PrintMatrix(self.transitionMatrix)

        for s in range(self.numReachableStates):
            print "ComputeTotalExpectedRewardWithAlgo2_SPIN2009(): self.selfPreXPost[%d] = %s." % (s, str(self.selfPreXPost[s]))
            sys.stdout.flush()

            if ((s == s0) or (s == sB)):
                continue

            for (s1, s2) in self.selfPreXPost[s]:
                pe = self.transitionMatrix[s1][s] * (1 / (1 - self.transitionMatrix[s][s])) * self.transitionMatrix[s][s2]
                re = self.rewardTransition[s1][s] + self.rewardTransition[s][s2] + self.rewardState[s] + self.transitionMatrix[s][s] / (1 - self.transitionMatrix[s][s]) * (self.rewardTransition[s][s] + self.rewardState[s])

                pe = sympy.simplify(pe) #, SYMN
                re = sympy.simplify(re) #, SYMN

                print "pe =", pe
                print "re =", re

                #re = self.reward[s1][s] + self.reward[s][s2] + self.rewardState[s] + self.transitionMatrix[s][s] / (1 - self.transitionMatrix[s][s]) * (self.reward[s][s] + self.rewardState[s])
                self.rewardTransition[s1][s2] = (pe * re + self.transitionMatrix[s1][s2] * self.rewardTransition[s1][s2]) / (pe + self.transitionMatrix[s1][s2])
                self.transitionMatrix[s1][s2] += pe

                self.rewardTransition[s1][s2] = sympy.simplify(self.rewardTransition[s1][s2]) #, SYMN
                self.transitionMatrix[s1][s2] = sympy.simplify(self.transitionMatrix[s1][s2]) #, SYMN

            print "We now remove %d from pre, post, self.selfPreXPost and put together the edges, etc." % s
            """
            #We remove s from the self.selfPreXPost[] lists. Then, we connect the pre and post vertices of s directly (result of the removal of s)
            for myIndex in range(self.numReachableStates):
                if (myIndex == s):
                    continue

                i = 0
                while (i < len(self.selfPreXPost[myIndex])):
                #for i in range(len(self.selfPreXPost[myIndex])):
                    if (self.selfPreXPost[myIndex][i][0] == s):
                        #del self.selfPreXPost[myIndex][i]
                        self.selfPreXPost[myIndex].pop(i)
                        #self.selfPreXPost[myIndex].pop(i)
                        try:
                            self.pre[self.selfPreXPost[myIndex][i][1]].remove(s)
                        except:
                            pass
                    elif (self.selfPreXPost[myIndex][i][1] == s):
                        self.selfPreXPost[myIndex].pop(i)
                        #self.selfPreXPost[myIndex].pop(i)
                        try:
                            self.post[self.selfPreXPost[myIndex][i][0]].remove(s)
                        except:
                            pass
                    else:
                        i += 1

                print "After remove s=%d, self.selfPreXPost[%d] = %s." % (s, myIndex, str(self.selfPreXPost[myIndex]))
            """

            #We connect the pre and post vertices of s directly (result of the removal of s). Then, we recompute the self.selfPreXPost[] lists.
            for i in range(len(self.pre[s])):
                for j in range(len(self.post[s])):
                    if (i < len(self.pre[s])): #The remove(s) below might violate these constraints.
                        sPre = self.pre[s][i]

                        if (j < len(self.post[s])): #The remove(s) below might violate these constraints.
                            sPost = self.post[s][j]

                            if (sPre != sPost): #For ComputeTotalExpectedRewardWithAlgo2_SPIN2009 to work we don't consider as pre and post the state itself, if it has a self-loop.
                                #We append sPost to self.post[sPre], if it doesn't exist already
                                try:
                                    res = self.post[sPre].index(sPost)
                                except:
                                    #pass
                                    self.post[sPre].append(sPost)

                        try:
                            res = self.post[sPre].remove(s)
                        except:
                            pass

                    if (j < len(self.post[s])): #The remove(s) below might violate these constraints.
                        sPost = self.post[s][j]

                        if (i < len(self.pre[s])): #The remove(s) below might violate these constraints.
                            sPre = self.pre[s][i]

                            if (sPre != sPost): #For ComputeTotalExpectedRewardWithAlgo2_SPIN2009 to work we don't consider as pre and post the state itself, if it has a self-loop.
                                #We append sPre to self.post[sPost], if it doesn't exist already
                                try:
                                    res = self.pre[sPost].index(sPre)
                                except:
                                    #pass
                                    self.pre[sPost].append(sPre)

                        try:
                            res = self.pre[sPost].remove(s)
                        except:
                            pass

            self.ComputeSelfPreXPost()

            self.PrintMatrix(self.transitionMatrix)
            self.PrintMatrix(self.rewardTransition)

        res = self.transitionMatrix[s0][sB] / (1 - self.transitionMatrix[s0][s0]) * (self.rewardState[s0] + self.rewardTransition[s0][sB])
        res += self.transitionMatrix[s0][s0] * self.transitionMatrix[s0][sB] / ( (1 - self.transitionMatrix[s0][s0]) * (1 - self.transitionMatrix[s0][s0]) ) * (self.rewardState[s0] + self.rewardTransition[s0][s0]) 
        #print "ComputeTotalExpectedRewardWithAlgo2_SPIN2009(): res before simplify =", res
        res = sympy.simplify(res) #simplify(res, SYMM)
        #sympy.simplify(res, SYMN)

        print "ComputeTotalExpectedRewardWithAlgo2_SPIN2009(): res =", res
        sys.stdout.flush()

        return res



###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################




mrm = None
SYMN = sympy.Symbol("N")
SYMM = sympy.Symbol("M") #We can substitute M with N - i.
SYMi = sympy.Symbol("i")


def Case1():
    #A CFG with if-then-else branches with probabilities 0.1 and 0.9.
    mrm = MDP(numStateVars = 1, stateVarsCardinality = [4], numActions = 1, initState = [0]) #This is an MRM

    mrm.transitionMatrix[0][0][1] = 0.1
    mrm.transitionMatrix[0][0][2] = 0.9
    mrm.transitionMatrix[0][1][3] = 1.0
    mrm.transitionMatrix[0][2][3] = 1.0
    mrm.transitionMatrix[0][3][3] = 1.0

    mrm.rewardState = [0.0, 1.0, SYMN, 0.0]


def Case2():
    #CFG of a loop: i = 0..N-1 . - putting symbolic value on transition probability ---> the average-time is a big 
    #    "series" which hardly? converges
    mrm = MRM(numStateVars = 1, stateVarsCardinality = [5], numActions = 1, initState = [0]) #This is an MRM

    mrm.transitionMatrix[0][0][1] = 1.0
    mrm.transitionMatrix[0][1][2] = 1.0
    mrm.transitionMatrix[0][2][3] = 1.0
    mrm.transitionMatrix[0][3][1] = (SYMN - 1) / SYMN
    mrm.transitionMatrix[0][3][4] = 1 / SYMN
    mrm.transitionMatrix[0][4][4] = 1.0

    mrm.rewardState = [0.0, 0.0, 1.0, 0.0, 0.0]


def Case3():
    #CFG of a loop: i = 0..N-1 without backedge and with multiplied costs of inner blocks with N, in order 
    #to avoid blowup of symbolic computation.
    mrm = MRM(numStateVars = 1, stateVarsCardinality = [5], numActions = 1, initState = [0]) #This is an MRM

    mrm.transitionMatrix[0][0][1] = 1.0
    mrm.transitionMatrix[0][1][2] = 1.0
    mrm.transitionMatrix[0][2][3] = 1.0
    mrm.transitionMatrix[0][3][4] = 1.0
    mrm.transitionMatrix[0][4][4] = 1.0

    c0 = 1
    c1 = 1
    c2 = 1
    c3 = 1
    c4 = 0
    mrm.rewardState = [c0, c1 * SYMN, c2 * SYMN, c3 * SYMN, c4]


def Case4():
    #CFG of insertion/selection? sort: 2 loops, one i = 0..N-1, the other j = i..N-1.
    #  for (i = 0; i < N; i++)
    #      for (j = i; j < N; j++)
    #          if a[j] < min ...
    mrm = MRM(numStateVars = 1, stateVarsCardinality = [7], numActions = 1, initState = [0]) #This is an MRM

    mrm.transitionMatrix[0][0][1] = 1.0
    mrm.transitionMatrix[0][1][2] = 1.0
    mrm.transitionMatrix[0][2][3] = 1.0
    mrm.transitionMatrix[0][3][4] = 1.0
    mrm.transitionMatrix[0][4][3] = (SYMN - SYMi - 1) / (SYMN - SYMi)
    mrm.transitionMatrix[0][4][5] = 1 / (SYMN - SYMi)
    mrm.transitionMatrix[0][5][2] = (SYMN - 1) / SYMN
    mrm.transitionMatrix[0][5][6] = 1 / SYMN
    mrm.transitionMatrix[0][6][6] = 1.0

    mrm.rewardState = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]


def Case5():
    #????
    #CFG of insertion/selection? sort: 2 loops, one i = 0..N-1, the other j = i..N-1.
    #  for (i = 0; i < N; i++)
    #      for (j = i; j < N; j++)
    #          if a[j] < min ...

    #Note that SYMi is actually not an input var, but the induction var of the outer loop: i = 0..N-1
    #Avg-time is 1.0*N*(0.5*N - 0.5*i) + 1.0*N*(1.0*N - 1.0*i) + 0.8*N = 1.5 * N * (N - i) + 0.8 * N
    mrm = MRM(numStateVars = 1, stateVarsCardinality = [7], numActions = 1, initState = [0]) #This is an MRM

    mrm.transitionMatrix[0][0][1] = 1.0
    mrm.transitionMatrix[0][1][2] = 1.0
    mrm.transitionMatrix[0][2][3] = 1.0
    mrm.transitionMatrix[0][3][4] = 1.0
    mrm.transitionMatrix[0][4][5] = 1.0
    mrm.transitionMatrix[0][5][6] = 1.0
    mrm.transitionMatrix[0][6][6] = 1.0

    mrm.rewardState = [0.0, 0.0, 0.3 * SYMN, 1.0 * (SYMN - SYMi) * SYMN, 0.5 * (SYMN - SYMi) * SYMN, 0.5 * SYMN, 0.0]


def Case6():
    global mrm

    #Note that the loop is modelled with loop exit condition checked at the footer (which is less realistic).
    #CFG of 2 loops
    #  for (i = 0; i < N; i++)
    #      for (j = 0; j < N; j++)
    #          if a[j] < min ...

    mrm = MRM(numStateVars = 1, stateVarsCardinality = [7], initState = [0]) #This is an MRM

    #Initially transition probabilities:
    mrm.transitionMatrix[0][1] = 1.0
    mrm.transitionMatrix[1][2] = 1.0
    mrm.transitionMatrix[2][3] = 1.0
    mrm.transitionMatrix[3][4] = 1.0
    mrm.transitionMatrix[4][3] = (SYMM - 1) / SYMM
    mrm.transitionMatrix[4][5] = 1 / SYMM
    mrm.transitionMatrix[5][2] = (SYMN - 1) / SYMN
    mrm.transitionMatrix[5][6] = 1 / SYMN
    mrm.transitionMatrix[6][6] = 1.0

    #mrm.PrintTransitionMatrix()

    mrm.rewardState = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]


def Case7():
    global mrm

    mrm = MRM(numStateVars = 1, stateVarsCardinality = [2], initState = [0])

    #Initially transition probabilities:
    mrm.transitionMatrix[0][0] = (SYMM - 1) / SYMM
    mrm.transitionMatrix[0][1] = 1 / SYMM
    #mrm.transitionMatrix[1][1] = 1.0

    mrm.rewardState = [1.0, 0.0]


def Case8():
    """
    For this case, the right result is: M
        returned by mrm.ComputeTotalExpectedRewardWithAlgo2_SPIN2009().
    However, mrm.ComputeTotalExpectedRewardWithValueIteration() returns:
        (9.0*M**8 - 36.0*M**7 + 84.0*M**6 - 126.0*M**5 + 126.0*M**4 - 84.0*M**3 + 36.0*M**2 - 9.0*M + 1.0)/M**8, which is close to O(M).
    """
    global mrm

    mrm = MRM(numStateVars = 1, stateVarsCardinality = [3], initState = [0])

    #Initially transition probabilities:
    mrm.transitionMatrix[0][1] = 1
    mrm.transitionMatrix[1][1] = (SYMM - 1) / SYMM
    mrm.transitionMatrix[1][2] = 1 / SYMM
    #mrm.transitionMatrix[2][2] = 1.0

    mrm.rewardState = [0.0, 1.0, 0.0]


def Case_InsertionSort():
    global mrm

    #Note that the loop is modelled with loop exit condition checked at the footer (which is less realistic).
    #CFG of insertion/selection? sort: 2 loops, one i = 0..N-1, the other j = i..N-1.
    #  for (i = 0; i < N; i++)
    #      for (j = i; j < N; j++)
    #          if a[j] < min ...

    mrm = MRM(numStateVars = 1, stateVarsCardinality = [7], initState = [0]) #This is an MRM

    #Initially transition probabilities:
    mrm.transitionMatrix[0][1] = 1.0
    mrm.transitionMatrix[1][2] = 1.0
    mrm.transitionMatrix[2][3] = 1.0
    mrm.transitionMatrix[3][4] = 1.0
    mrm.transitionMatrix[4][3] = (SYMN - SYMi - 1) / (SYMN - SYMi)
    mrm.transitionMatrix[4][5] = 1 / (SYMN - SYMi)
    mrm.transitionMatrix[5][2] = (SYMN - 1) / SYMN
    mrm.transitionMatrix[5][6] = 1 / SYMN
    mrm.transitionMatrix[6][6] = 1.0

    #mrm.PrintTransitionMatrix()

    mrm.rewardState = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]


###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################
###############################################################################

def Main():
    #f = lambda N: sympy.Sum(i, (i, 1, N))
    #print sympy.solve( sympy.sum(n, (n, 1, N)), N )
    #print sympy.solve( sympy.Sum(SYMi, (SYMi, 1, SYMN)), SYMN ) #ValueError: Could not parse expression Sum(i, (i, 1, N))

    #We need to compute sum of symbolic vars
    if (False):
    #if (True):
        n = sympy.var("n")

        print sympy.Sum(1/n**2, (n, 1, sympy.oo))
        print sympy.Sum(1/n**2, (n, 1, sympy.oo)).evalf()

        #print sympy.Sum(SYMi, (SYMi, 0, SYMN - 1)).evalf()
        print sympy.Sum(SYMi, (SYMi, 0, SYMN - 1)).doit()

        print sympy.summation(SYMi, (SYMi, 0, SYMN - 1))

        print sympy.nsimplify(sympy.Sum(1/n**2, (n, 1, sympy.oo)), [sympy.pi])

        #print sympy.nsimplify(sympy.Sum(SYMN - SYMi, (SYMi, 0, SYMN - 1)), [n]) #ValueError: constants must be real-valued
        print sympy.simplify(sympy.Sum(SYMN - SYMi, (SYMi, 0, SYMN - 1)))

        f = lambda N: sympy.Sum(n, (n, 1, N))

        x = sympy.var("x")
        print f(x)

        N = sympy.var("N")

        print sympy.solve(f(x)) #ValueError: Could not parse expression Sum(n, (n, 1, x))

        print sympy.solve( sympy.Sum(n, (n, 1, N)), N )

    #return

    #Case6()
    #Case7()
    #Case8()
    Case_InsertionSort()

    mrm.PrintMatrix(mrm.transitionMatrix)

    res = mrm.ComputeTotalExpectedRewardWithAlgo2_SPIN2009(0, mrm.numReachableStates - 1)
    #print dir(res)

    print "res =", res
    print "res.has(SYMN - SYMi) =", res.has(SYMN - SYMi)
    print "res.has(2*SYMN - SYMi) =", res.has(2*SYMN - SYMi)
    #!!!!We need to check if N is from outer loop and only then do sum(res / N). Generalize for more loops, etc.
    if (res.has(SYMi)):
    #if (True):
        print "Sum_{i=0}^{N-1}(res / N) =", sympy.Sum(res / SYMN, (SYMi, 0, SYMN - 1)).doit()
    """
    Other methods for res:
        __contains__
        is_Symbol
        match
        is_NumberSymbol
        free_symbols
        flatten
        factor
        apart
    """

    if (False):
    #if (True):
        mrm.ComputeTotalExpectedRewardWithValueIteration()


if __name__ == "__main__":
    if (False):
    #if (True):
        i = sympy.Symbol("i", real = True)
        N = sympy.Symbol("N", real = True)
        M = sympy.Symbol("M", real = True)

        print sympy.simplify(M*(1 - (M - 1)/M)**2)
        print sympy.simplify(1.0*(M - 1)/(M*(1 - (M - 1)/M)**2))
        print sympy.simplify(1.0/(1 - (M - 1)/M) + 1.0*(M - 1)/(M*(1 - (M - 1)/M)**2))
        #"""
        #simplify an expression - See file:///Z:/1PhD/ReVival/sympy/Docs/sympy-0.7.1-docs-html/modules/simplify/simplify.html#sympy.simplify.simplify.combsimp: "Collect additive terms with respect to a list of symbols up to powers with rational exponents."
        #print sympy.simplify
        print sympy.simplify(1.0*N*(0.5*N - 0.5*i) + 1.0*N*(1.0*N - 1.0*i) + 0.8*N, N)
        print sympy.collect(1.0*N*(0.5*N - 0.5*i) + 1.0*N*(1.0*N - 1.0*i) + 0.8*N, N)
        print "Am here"
        #"""

        print sympy.summation(i, (i, 1, N))
        print sympy.summation(N - i, (i, 0, N - 1))

    Main()
